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EIGENSYSTEM ANALYSIS OF CLASSICAL RELAXATION TECHNIQUES 
WITH APPLICATIONS TO MULTIGRID ANALYSIS 


Harvard Lomax, Catherine M. Maksymiuk 
Ames Research Center 


SUMMARY 

Classical relaxation techniques are related to numerical methods for solution of 
ordinary differential equations. Eigensystems for Point-Jacobi, Gauss-Seidel, and 
SOR methods are presented. Solution techniques such as eigenvector annihilation, 
eigensystem mixing, and multigrid methods are examined with regard to the eigen- 
structure. 


1. MATRIX FORMS OF FINITE DIFFERENCE SCHEMES 
1.1 Banded Matrices 

The symbol B[M : arguments) is used to represent a matrix of order M, the 
elements of which are all zero except for those along diagonals close to the central 
diagonal. The number of arguments is always odd and the central argument refers 
to entries in the central diagonal. The M is often omitted from the argument list. 
Thus 


B(b_2,b..i,bo,b\,b2) ( 1 - 1 . 1 ) 

represents a matrix with scalar entries that are constant along the central diagonal 
and the two diagonals just above and just below it. This is a pentadiagonal ma- 
trix. A form of banded matrix that is very common in numerical analysis is the 
tridiagonal system illustrated by 


1 


b i c 2 
ai 62 C3 

0-2 bz 


B M 


: a, 6, c j 


( 1 . 1 . 2 ) 


cm 

Ojw-l &mJ 

In this case the elements vary along the diagonals and are represented by vectors. 
The particular indexing of the elements, constant along columns, is used in so-called 
conservation-law formulations. 

Banded matrices can be used to represent finite difference schemes, in which 
case the nature of the boundary conditions is signified by the matrix formulation. 
Eq (1.1.2) would be used if the boundary conditions were Dirichlet. If periodic 
boundary conditions were desired, the appropriate form would be 

T b\ c 2 a-M 

0\ 62 c 3 


B p (m : a,b,cj = 


(1.1.3) 


An important set of banded matrices with constant diagonal entries is repre- 
sented by the following expression 


B r (a, b f c) 


(1.1.4) 


2 



This particular matrix is referred to as a tridiagonal circulant matrix. Circulant 
matrices need not be tridiagonal; they can be completely dense. Notice that a 
circulant matrix is a special form of a periodic matrix. 


The notation I is used for the identity matrix and the notation D(b), is used 
for a diagonal matrix with constant elements. Notice that 


B{b) = D{b) and B{ 1) = £>(1) = / 


( 1 . 1 . 5 ) 


1.2 Difference Schemes as Banded Matrices 


A difference scheme is usually written as a point operator, 
central-difference scheme for a second derivative is given by 


The three-point 


— ^^ 2 ( u .7 + 1 T u j- 1) 


( 1 . 2 . 1 ) 


where the index j refers to the location of the dependent variable in the equispaced 
x direction. This difference scheme can be expressed as a matrix operator 
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if the boundary conditions are speci- 
fied. If u is specified at a boundary, the 
(6c) (boundary condition) is referred to 
as Dirichlet, and if is specified at 

a boundary, the (be) is referred to as 
Neumann. Figure (1.2.1) shows three 
possibilities. In Fig. (1.2. la) a mesh 
with 4 interior points is shown. The 
value of u is specified at the two end 
points a and 6, so the (6c) are said 
to be Dirichlet. Figure (1.2.1b) illus- 
trates a situation representing a Neu- 
mann (6c) for the right side and a 
Dirichlet (6c) for the left side. Periodic 
(6c) are illustrated in Fig. (1.2.1c). No- 
tice that the normalized length of the 
mesh is n in the top case, 7r/2 in the 
middle case, and 2n in the bottom case. 
These conventions fix a condition be- 
tween the space step size and the num- 
ber of points in a mesh which is conve- 
nient in later developments. 

The matrix representation of eq (1.2.1) 
for the three kinds of (6c) illustrated 
can now be written. For the com- 
plete Dirichlet problem illustrated in 
Fig. (1.2.1a) 


X = 0 Ax = 7r/(M + 1) X = 7T 


h 


4 

M 


a. Dirichlet (6c) on both sides 


7T 


x = 0 Ax = 7 t/( 2M + l) x = — 

2 


0 

1 
1 


o 

2 


H 

o 

3 


o 

4 6 

M 


b. Neumann (6c) on right 


x = 0 Ax = 27 r/M x — 2 tx 


o 

4 


6 

1 


o 

2 


o 

3 


o 

4 


6 

1 


c. Periodic (6c) 


Figure 1.2.1 - Boundary conditions. 


6xxU = X^[b( 1,-2,1)u+ (6c)] 
(6c) [u a , 0, • • • , 0, uj,] 


For the mixed Neumann-Dirichlet problem in Fig. (1.2.1b) 
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— 


6 = 


(be) = 


1 


Ax 2 

[ 2 , 2 , 2, 1 


B^l, — 6, l^u + ^6c^ 

T 
T 


( du\ 

a,°,0,Az(^ — J^j 


(1.2.3) 


It should be noted that for the Neumann (6c) the value of is assumed to given at 
the point M + §Az. Futhermore this particular scheme is only first order accurate, 

although the multiplying constant is small. Finally, for periodic (6c) shown in 
Fig. (1.2.1c) we have the relation 


rU - 


Ax 2 


B P ( l,-2,l)u 


(1.2.4) 


Notice that the periodic case has no (6c) and the expression is homogeneous. 
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2. FORMULATION OF THE MODEL PROBLEM 


2.1 The Basie Equation and its Solution 

Is is assumed that some time dependent partial differential equation represents 
some valid fluid flow. It is further assumed that this equation has proper boundary 
conditions and that appropriate differencing schemes have been applied to approxi- 
mate the space derivatives and the boundary conditions. The result is representable 
in matrix form as 


au — , 

— + A b u - f h = 0 (2.1.1) 

where u represents the dependent variables, f b holds the boundary conditions and 
the forcing function (if there is one), and At, is usually nonlinear (i.e., depends on 

u). Now it is assumed that a solution to this equation exists for which u is time 

invariant. We refer to this as a steady state solution of eq (2.1.1). Such a solution 
satisfies the equations 


A b u-f b = 0 (2.1.2a) 

or u = A^ l f b (2.1.2b) 

where A^ 1 is assumed to exist. None of the questions of existence and well-posed 
boundary conditions are considered in this paper. Finally, it is assumed that the 
solution given by eq (2.1.2) is going to be found by an iterative process, and we 
wish to study methods for carrying out these iterations. 

The above was written to treat the general case. It is instructive in formulating 
the concepts to consider the special case given by the linear diffusion equation in 
one dimension 


du 

dt 


d 2 u 
dx 2 


p{x) 


This has the steady state solution 


(2.1.3) 


d 2 u 
dx 2 


p{x) 


(2.1.4) 


which is the one dimensional form of the Poisson equation. Introducing the three- 
point central differencing scheme for the second derivative, eq (1.2.1), we find 
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du 


1 


B( 1, -2,1 )u + (fee) - p 


dt Ax 2 

where (fee) contains the boundary conditions. In this case 


(2.1.5) 


A* = ^1,-2,!) 

h = P~ (he) 


( 2 . 1 . 6 ) 


We work extensively with this model form. 

2.2 Preconditioning the Basic Matrix 

The iteration process mentioned above is referred to as a relaxation procedure. 
It is standard practice in applying relaxation procedures to precondition the basic 
equation. This preconditioning has the effect of multiplying eq (2.1.2a) from the left 
by some nonsingular matrix. In the simplest possible case the conditioning matrix 
is a diagonal matrix composed of a constant D(b), see eq (1.1.5). If we designate 
the conditioning matrix by C, the problem becomes one of solving for u in 


CA b u - Cf b = 0 


( 2 . 2 . 1 ) 


Notice that the solution of eq (2.2.1) is 

u = \CA b )-'Cf h = A b 1 C~ 1 Cf b = V/t (2-2.2) 

which is identical to the steady-state solution of eq (2.1.1), provided C -1 exists. 

In the following we will see that our approach to the iterative solution of eq 
(2.2.1) depends crucially on the eigenvalue and eigenvector structure of the matrix 
\C A b \, and, equally important, does not depend at all on the eigensystem of the 
basic matrix A b . A simple example illustrates the point. Consider the first-order 
backward difference scheme (If properly applied this can also be referred to as an 
“upwind” scheme). 


(6x u ) = ^( U J ~ u J-i) (2.2.3) 

If u is specified on the boundary at the left ( u a given in Fig. 1.2.1a), the difference 
scheme forms the matrix operation 
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S x u = £?(-!, l,0)u + 


(2.2.4) 


With this operator, the simple relation 

g(x\)dx] (2.2.5) 

can be approximated by 

2-B(-I,l,0);+^(£) =»' (2.2.6) 

where (bc^j = |-u a ,0, ,0] T 

Eq (2.2.6) takes the form of eq (2.1.2a) if 

M, ~ ^B(- 1,1,0) and -^ = 9 - Xx( 6c ) (2.2.7) 

The eigensystem of At, is fully defective. 1 However, let C = — Ax- 5(0, 1, -1). The 
resulting product matrix \CA}\, see eq (2.2.9), has a complete set of eigenvectors 
and the eigenvalues are all real and negative. The relaxation of eq (2.2.1) with this 
choice of C and At, can proceed along well defined classical lines in spite of the fact 
that the original basic matrix had no resemblance to the classical form. 

Finding the solution to eq (2.2.6) is a trivial matter, but that is not our object. 
Our object is to use eq (2.1.2) and (2.2.7) to illustrate a general iterative solution 
process. This process consists of first, choosing C so that [CMf,] is “close” to a model 
matrix given in the next section, and second, making a thorough analysis of the 
various iteration techniques available for the model matrices. There are well-known 
techniques for accelerating relaxation schemes if the eigenvalues of \C Ab\ are all real 
and of the same sign. This paper is limited to a study of only these techniques. 
Such a limitation imposes the requirement: 

The conditioning matrix C is chosen such that the eigen- 
values of [CMt] are all real and negative, that is, such 
that \CAt\ is negative definite. (2.2.8) 

*See Section 4.7 


= </(*) 


“= f 

JO 
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It is also well known that a choice of C which ensures this condition is the neg- 
ative transpose of Ab . In the above example given by eq (2.2.7), B T ( 1,1,0) = 
B( 0. 1,-1) and 


-5(0,1, -1)5(-1, 1,0) 


- 2 1 
1-2 1 
1 - 2 1 
1-2 1 
1 - 1 


(2.2.9) 


which is the same matrix that arises for a second-derivative approximation with 
mixed Dirichlet-Neumann boundary condition, see eq (1.2.3) and (5.6.2). It is easy 
to verify that if, instead of C — —B{ 0, 1,-1), C is chosen to be 


C = 


1 - 1 
1 - 1 
1 - 1 
1 - 1 

11112 


( 2 . 2 . 10 ) 


the product C ■ B(- 1,1,0) is equal to 5(1, -2,1), a finite-difference matrix for a 
second derivative with Dirichlet conditions on both sides, see eq (5.6.1). 


Another interesting example arises from the study of first-order partial differ- 
ential equations when a central differencing (rather than upwind, as in eq (2.2.3)) 
scheme is used for the approximation of the derivative. However, the physics of 

this problem permits a Dirichlet on one end but allows no constraint on the 

other; so on one end (and we assume it is the right side) we use a first order upwind 
scheme. This leads to the matrix difference equation 


6 x u — 


1 

2Ax 



0 1 


Ua 



- 1 0 1 


0 


< 

- 1 0 1 

u + 

0 



- 1 0 1 


0 


< 

- 2 2 


0 



( 2 . 2 . 11 ) 


The matrix in eq (2.2.11) can first be conditioned so that the modulus of each ele- 
ment is 1, and then further conditioned with multiplication by the negative trans- 
pose. The result is 
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a 2 = -AjAl = 


0 1 


0 1 

- 1 0 1 


- 1 0 1 

- 1 0 1 


- 1 0 1 

- 1 0 1 


- 1 0 1 

- 1 - 1 


- 1 1 


-10 1 

0-2 0 

1 0-2 
1 0 

1 


1 

0 1 

2 ] 

1 - 2 


( 2 . 2 . 12 ) 


The matrix on the right side of eq (2.2.12) does not look familiar, but if we define 
a permutation matrix P and carry out the process P T \~Aj A \}P we find 


'o 

1 

0 

0 

o" 


' - 1 

0 

1 




'o 

0 

0 

0 

r 

0 

0 

0 

1 

0 


0 

- 2 

0 

1 



1 

0 

0 

0 

0 

0 

0 

0 

0 

1 


1 

0 

- 2 

0 

1 


0 

0 

0 

1 

0 

0 

0 

1 

0 

0 



1 

0 

- 2 

1 


0 

1 

0 

0 

0 

1 

0 

0 

0 

0 




1 

1 

- 2 


0 

0 

1 

0 

0 


2 1 

1-2 1 
1-2 1 
1-2 1 
1 - 1 


(2.2.13) 


which leads to exactly the same matrix as that derived from upwind differencing 
presented in eq (2.2.9). 


The importance of these concepts is much more evident when they are used 
to precondition the Cauchy-Riemann and Euler equations. In these cases even 
when the basic matrix At, has nearly imaginary eigenvalues, as it does if central 
difference schemes are used for the first, derivatives, the conditioned matrix [-T^\4f,j 
is nevertheless negative definite and the study of the model equations in the next 
section is pertinent to its solution. 
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2.3 The Model Equations 


Preconditioning processes such as those described in the last section allow us 
to prepare our algebraic equations in advance so that certain eigenstructures are 
guaranteed. In the remaining part of this paper we wish to thoroughly investigate 
some simple equations that model these structures. There is absolutely nothing 
new in these model equations; they (especially the two-dimensional case) are used 
in all the standard texts treating the subject of relaxation. However, the approach 
to their solution, i.e., as a subset of the theory of ordinary differential equations, 
is novel in some places, and the picture of the subject from this point of view is 
helpful in unifying a large number of relaxation methods that have been proposed. 

Consider the preconditioned iterative equation having the form 

^ + A* - / - 0 (2.3.1) 

where A is negative definite and the symbol for the dependent variable has now been 
changed to <f> signifying that the physics being modeled is no longer time accurate. 
However, notice that a steady state of eq (2.3.1) is guaranteed to exist (because it 
is negative definite) and that steady state solution is A -1 /, which is identical to eq 
(2.1.2). In the notation of eq (2.1.2) and (2.2.1) 

A = CA h and f = Cf h (2.3.2) 

For the model equations, in the one-dimensional case A has the form 

A = fi(l,6,l) 

b= [ — 2 , - 2 , , -2,s} T 

s = —2 or — 1 (2.3.3) 

For the two -dimensional case, A has the form 

A B^oiyl , B{ol x , 4, ot x ) , Oy/], 

ot x A ct y = 2 ; a x < 2 (2.3.4) 

A tremendous amount of insight to the basic features of relaxation is gained 
by an appropriate study of the one-dimensional case, and much of the remaining 
material is devoted to this case. We attempt to do this in such a way, however, that 
it is directly applicable to two- and three-dimensional problems. 
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3. THE DELTA FORM OF AN ITERATION SCHEME 


3.1 Basic Theory 

We assume that our basic difference matrix has been preconditioned to form 


A4> - f = 0 (3.1.1) 

and thatA is negative definite. We choose some iterative scheme to find the solution, 
A -1 /* and designate the iteration count by the subscript n or the superscript (n). 
The converged solution is designated <f> ^ so that 

= (3.1.2) 

The manner in which convergence is measured in actual practice is not considered 
here. (It is generally related to the magnitude of a residual, e.g. £ \ A<j>- f\, summed 
at each point in the mesh.) Instead we assume that after some step N , the solutions 
,4>n + \,- • -,<j>N+k are all “close enough” that each could be considered the final 
answer. That is to say, for our purposes, one can write 4>s = <i>N + 1 = • ■ • = 4>N-rk- 
Suppose F(z) is some function having the property .F(O) = 0. Then the general 
expression for the “delta form” used to relax eq (3.1.1) can be written 

- $ n+l ) = At. - h { 1 : (*-»■*> 

If F is a linear operator, one can easily prove by the theory of finite difference 
equations that the particular solution of eq (3.1.3) is our desired solution A -1 /, 

if f is independent of n. It is only necessary then to formulate F in such a way 
that the complementary solution of the difference equation (3.1.3) goes to zero as 

n — ► oo. 


3.2 Examples 


Consider the model equation 

-6(1, -2, 1)^ -/ = 0 (3.2.1) 


One example of a delta form, written as a point operator, is 
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«(^" +1 > - <-") + f>(*<" +1) - 2*<"> + = *<*’, - 2*‘"> + *<?, - f, 

(3.2.2) 

Notice that the left side sums to zero if <f>^ converges. If we set b — 1 and 
a ~ 1/(2 h), eq (3.2.2) reduces to 


,(n+l) A^-i) 


= 2 fc («‘” ) 1 - ^" +1) - 11 + - />) 


(3.2.3) 


which is a representation of the D u Fort- Fran kel method used to solve the diffusion 
equation, see e.g. ref. 1, p. 60. Notice also that eq (3.2.2) can be interpreted in the 
form of an ordinary differential equation 


+ 4t= e(i ’ 


(3.2.4) 


if the iteration index is thought of as a “time” displacement. 

A second example of a delta form written as a point operator for eq (3.2.1) is 

- *< n > - + *£.’,) + ^ n+i) - ^■ , ) = - 2 ^ - 1 , 

(3.2.5) 

Again we see that the left side sums to zero at convergence. This time, however, 
the expression can be interpreted in the form of a partial differential equation 


d 2 <t> , u d<j> _ d 2 4 > A 

+ 6 at “ a? _ /(l) 


(3.2.6) 


This approach to relaxation has been used by Garabedian, see e.g. ref. 1, p. 125. 
3.3 The Ordinary Differential Equation Formulation 

The particular type of delta form considered in the remaining part of this 
paper is one that leads to a differential interpretation composed of a set of coupled, 
first-order, ordinary, differential equations. In difference notation it is expressed as 


H 4> n + x <p n = A<j> n - f 
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where H is some nonsingular matrix which is independent of n for “stationary” 
methods and is a function of n for “nonstationary” ones. In differential notation it 
can be written 


H -£ =A +n-f* ( 3 ‘ 32) 

Since the basic matrix Ai, has already been conditioned by the matrix C to 
produce A = C Ai and / = C f h , the secondary conditioning matrix H may seem 
superfluous. We find it convenient, however, to separate the relaxation procedure 
into two distinct parts. First, we precondition the basic equation to put it in a 
model lorm. Second, we lurther condition the model form to generate optimum 
algorithms. For this reason both C and H play a useful role. 
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4. THE ANALYTICAL AND NUMERICAL SOLUTION OF FIRST-ORDER 
ORDINARY DIFFERENTIAL EQUATIONS WITH CONSTANT 
COEFFICIENTS AND A CONSTANT FORCING FUNCTION 


4.1 Introduction 

One purpose of this report is to present the subject of relaxation as a subset 
of the theory of the numerical solution of ODE. That concept in itself is certainly 
not new. However, a systematic treatment of such an approach has not, to the 
author’s knowledge, been published, and it leads to some interesting and useful 
interpretations. In order to make the discussion clear, a review of the theory for 
the analytical and numerical solution of ODE is given. 

4.2 The Analytical Solution of ODE 

We consider next the analytical solution of a set of coupled first-order, ordinary 
differential equations given by 


du 

dt 


[A]« - 7 


(4.2.1) 


where both [A] and / are independent of u and t. We first consider only those 
cases for which the eigenvectors of [A] are linearly independent. Introduce the left 
and right eigenvector matrices X~ l and A' such that 

Y" 1 AA r = £>(a) = A (4.2.2) 

where A is a diagonal matrix having the eigenvalues of [A] as the entries. When 

we speak of a single eigenvector of [A] , say x m , we are referring to a column in X 
corresponding to a A m in A. Linear independence of the eigenvectors means that 
a ■ x m + b • x n ± Xk for any values of (complex) a and b, and where m ± n ± k. To 
solve eq (4.2.1) we multiply each term from the left by X~ l and insert the identity 
matrix X~ l X = / between [A] and u. There results 

X" 1 ^ =X~ 1 AXX~ 1 u- X- 1 } (4.2.3) 

C hv 


which reduces to 


PRECEDING PAGE BLANK NOT FltMTO 
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Define 




(4.2.4) 


w = X _1 u ; g = X -1 / (4.2.5) 

and the solution can be written 

^ - D(A)w - g (4.2.6) 

Eqs (4.2.1) and (4.2.6) are expressing the same equality but in different algebraic 
forms. The key difference is that eq (4.2.6) is completely uncoupled. It can be 
written line by line as a set of first-order differential equations (defining w' as ^j>) 


w\ = Aitni - gi 

w' 2 = A 2 U> 2 - 02 

W rn ~ ^m w m ~ 9m | 


(4.2.7) 


each of which can be solved separately. The resulting solutions can then be recou- 
pled, using the inverse of eqs (4.2.5), to form the solution to eq (4.2.1). When the 

forcing function, g = X -1 /, is not a function of f, the solution of the mth equation 
in (4.2.7) is 


w m — c m e T « 9m (4.2.8) 

Tracing the algebra backward, one can easily show that the solution to eq (4.2.1) 
can be written 


u = X(ce A * } + XA -'X' 1 / 


(4.2.9) 


where ( ) denotes a vector of the enclosed terms. Thus 
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(4.2.10) 


t t\ — f A,t „ „ 0 

^C6 J — C\€ 9^2^ 9 


A M tl T 


Another way to express the solution given by eq (4.2.9) is 


u = cie*‘ t xi + C2e* 2< X 2 + ■ ■ • + cm e^ M< xm + [A] / 

Transient solution Steady state 


(4.2.11) 


Here in classical terminology [A] 1 / designates the particular solution to the com- 
plete equation and the remaining terms are the complementary solution to the 
homogeneous equation, du/dt = [A]u. In fluid flow problems [A] /is often 
referred to as the steady state solution and the complementary solution is often 
referred to as the transient solution. 

Given a linear system, the two possible formulations of the same problem 
discussed above and given by eqs (4.2.1) and (4.2.7) play an important role in 
our following development. For that reason we give them special names so we can 
quickly extract the concept in a given situation. This is summarized in the following: 


J t = [A]u - f{t) (4.2.12a) 

is referred to as an equation in real space, 


(4.2.12b) 

is referred to as an equation in eigen space (also referred to as wave space ) 

4.3 The Isolation Theorem and the Representative Equation 

Let us summarize the results of the previous section in order to simplify the 
discussion in the next section which treats the numerical solution of ODE. The 
fact that eqs (4.2.1) and (4.2.7) express the same equality is the basis for the 
Isolation Theorem which is stated next. 

Given a set of first order ordinary differential equations such that 


dw , . , - . 

— = [Aju; -g{t) 
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o The coupled equations are linear 
with time-constant coefficients. 

o The eigensystem of the [A] matrix 
is not defective. 


one can develop: 


The Isolation Theorem 1 


Applying any standard time- march method 
to each equation in a coupled set of 
equations having the constraints specified 
in (4.3.1) is mathematically identical to: 

o Uncoupling the set including the forcing terms, (4.3.2) 

c Individually integrating each equation in the 
uncoupled set, 

o Recoupling the group to form the final result. 

In the terminology of (4.2.12) this amounts to the observation that (under the 
conditions stated in (4.3.1)) using any time-march method on the coupled equations 
in real space is identical to using the same method on the single equation in eigen 
space. 

The theorem uses the terminology “mathematically identical to” which is rigor- 
ously correct. Unfortunately it is not strictly correct for the statement “numerically 
identical to” because of roundoff error. Numerical experiments with simple eigen- 
systems are easy to construct and quite informative in verifying the substance of 
the theorem. 

The “single equation in eigen space” , mentioned above, has the form shown by 
eq (4.2.7). We simplify this to 


(4.3.3) 

where A and a are (complex) constants. Notice that we only consider the case for 
a forcing function that, is independent of time, since that is sufficient when our 
application is in the field of relaxation. We refer to eq (4.3.3) as the re presentative 
equation. 


dw 

— = Aw + a 

dt 


1 Proof of this theorem is a simple exercise using the concepts 


outlined in section 


a. - + 
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4.4 The Numerical Solution of ODE 


Now we must choose some numerical procedure to integrate the representative 
eq (4.3.3), knowing from the isolation theorem that that method will advance each 
eigenvector in eq (4.2.11) according to its associated eigenvalue. These methods are 
referred to as time-marching methods and they convert the representative ODE into 
an ordinary difference equation (OAE), or a set of coupled OAE, depending upon 
the choice of method. These OAE can be solved analytically and their solution 
provides the result required to evaluate the method and compare it with other 
time-marching methods. 

First let us review the process for finding the solution of coupled first-order 
OAE that are linear and have constant (with respect to n) coefficients. These can 
be expressed in the form 


u n+ i = [C]u n -7 (4.4.1) 

following a procedure exactly the same as that used to express eq (4.2.1) in the 
form of eq (4.2.7), we can re-express eq (4.4.1) as 


(tt>n+l)l = 0\ {^n) 1 - 5\ 

(u>n+l)2 = 02{Wn)l ~ 5l 

(®n+l)m = Gm[ u ’n)m ~ Qm 


(4.4.2) 


where the a are the eigenvalues of [C] (which is assumed to be nondegenerate). 
This represents a set of uncoupled first-order OAE each of which has the form 

w n+ i=ow n ~g (4.4.3) 

This simple first order difference equation has the solution 

w n = c(o) n - — 1 — (4.4.4) 

1 — o 

which can be verified by substitution. 

Let us consider an example of how these results can be used to analyze a 
numerical time-marching method. We introduce some notation for the discrete 
variables, thus 
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(4.4.5) 


t = t n = nh ; h = At 
u n+k = u(t n + A:/i) 

u n+e = 

The numerical time-march method is described by a linear combination of the terms 
u n4 * and ti' n+<! where k,i = ± 0 , 1 , 2 ,... 

The simplest example of a numerical time-marching method is the explicit 
Euler scheme which can be written 

«n+ 1 - Un + hu' n 

Applying this to the coupled set of ODE represented by eq (4.2.1) gives 

u n + 1 ^77 T h A j u n — hf 

= [/ + h[A]}u n -hf 


(4.4.6) 


(4.4.7) 


Comparing this with eq (4.4.1) we see that for the explicit Euler method 

(4.4.8) 

First of all we notice that, since the identity matrix [/] commutes with any matrix 
and h is a scalar, the eigenvectors of [C] and [A] are identical. From this it follows 
that 


[C] = [7 + hk] 

f = hf 


c m — 1 T ^mh (4.4.9) 

Further, it can be shown (combining eqs (4.4.8) with eq (4.4.3) and recoupling the 
systems) that the particular solution of the OAE produced by the explicit Euler 
scheme is 


Wn -- C] ((7i) n Xi + C2(cr2) n a:2 T • • • + Cm(<7m)”^A/ + [A] f (4.4.10) 


where [A] l f is the exact solution of the ODE, and the eigenvectors x m are the 
eigenvectors of [A]. 
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4.5 The Analytic Solution of OAE 


For completeness we discuss briefly the analytic solution of linear difference 
equations with constant coefficients. For background, we recall that one classical 
approach to the study of linear ODE is carried out in terms of the operator D where 


D n 


d n 

dt n 


(4.5.1) 


The basic part of that solution process consists of replacing the time derivatives 
w ith the appropriate power of D and thereby deriving a characteristic polynomial, 
denoted P(D ) . This is followed by finding the roots, A m , of P(A) = 0 . The 
homogeneous part of the solution is then characterized by the expression 

u{t) = J2 c me X "' t (4.5.2) 

m 


where the c m are determined from the initial conditions. If u(t) is a vector, the 
terms on the right would be multiplied by the associated eigenvectors x m as in eq 
(4.2.1). 

A similar process can be carried out for linear OAE , except that the D opera- 
tor is replaced by the operator E, referred to as the displacement or shift operator , 
and defined in terms of the D operator by 

E = e hD = l + hD+ \h 2 D 2 + ••• (4.5.3) 

From this definition it should be clear that 

u n+ k = E k u n b n+a - E a b n u n + i = E?u n (4.5.4) 

where k and a can be integer, fractional, or even irrational numbers. Again the 
basic part of the solution process consists of finding a characteristic polynomial, 
this time designated P{E), and then finding the roots of P{o) — 0 . In this case for 
OAE the homogeneous part of the solution is characterized by 

= ]Tc m (<7 m ) n (4.5.5) 

m 

where o m are the roots of the characteristic polynomial and the c m are determined 
by the initial conditions. If u n is a vector, the terms on the right would be multiplied 
by the associated eignevectors as in eq (4.4.10). 
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As a simple example, consider the process for finding the solution to the equa- 
tion 


u n+2 + a l u n+\ + C|)U n = 6 n + a + 6” 

(4.5.6) 

In operational form this becomes 


P(E)u n = Q(E)b n 

(4.5.7) 

where 

■B(jE/) = E ^ 4" a\E -t Co 
Q{E) = E* + 1 

(4.5.8) 


We refer to Q(E) as the particular polynomial since it serves to determine the 
particular solution. The general solution of eq (4.5.7), which can be derived using 
Boole’s first rule (see ref. 2), is 


u n = c ™(°m) n + b n (4.5.9a) 

which reduces to 

U n = C rn{°m) n + ^7^ (4.5.9b) 

if the right side is independent of n (i.e. if b = 1). The general solution for the 
polynomials given by eq (4.5.8) is 


Un _ ^1 


a j -+• a j - 4a, 


+ c 2 


-ai - y'af - 4 go 


b a + 1 


_). ^ 1 J ' foTl 

6 2 + aj6 + do 

(4.5.10) 


4.6 The Analysis of Time-Marching Methods 

One can apply the analysis in Section 4.5 to investigate all manner of time- 
marching methods (e.g. Runge-Kutta, predictor corrector, and multistep methods) 
as they apply to the representative equation (4.3.3). Consider the first-order explicit 
Euler method for a scalar equation 

Un+i = u n + hu’ n (4.6.1) 

Applied to the representative equation it gives 

u n+i — (1 4 A h)u n 4 ha (4.6.2) 
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(4.6.3) 


and this produces eq (4.5.7) where 

P(E) = E - 1 - A h 
Q{E) - ah 

Using the results of eq (4.5.9b), we find 

u n = c(l + A h) n - ^ (4.6.4) 

Notice that the steady-state solution is the exact steady-state solution of the ODE 
and the o root approximates the Taylor series expansion of (e Xh ) through the order 
of the method. 

This o, A relation can be established for any linear time-march scheme such 
as Adams- Moulton and Runge Kutta, see e.g., ref. 3. It is fundamental to the ap- 
proach we present later because the convergence of the relaxation methods depends 
on |ct|, and the value of o is some function of the product A h. The exact nature of 
the functional dependence of o on A h is determined by the choice of the differencing 
scheme which fixes the characteristic polynomial P[E). If the time- march method 
produces only one o for each A, it is referred to as a one-root method. Many meth- 
ods generate more than one o for each A, in which case one of the o approximates 
e Xh and the others are referred to as spurious. The subject of spurious roots is 
outside the scope of this report. Some examples of the o, A relation are given in 
Table 4.6.1. 


TABLE 4.6.1 - THE o, A RELATION FOR SOME WELL-KNOWN 
NUMERICAL METHODS 


Method 


a, A relation 


Runge-Kutta of Order: 
First 
Second 

N th 

Leapfrog 

3-pt Adams-Bashforth 

Implicit Euler 
Trapezoidal method 


o — 

1 T A h 

a = 

l + \h+ \\ 2 h 2 

o — 

E 


°U2 = 

\h±y/\ + X 2 h 2 

<?1,2 = 

1 

2 

1 + ^ A hi yJ . 

a = 

1/(1 - A h) 

o — 

(1 

+ ^ A/i) / (l - 


l + A h+ |A 2 h 2 
\\h) 
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The importance of the cr, A relation is brought out very clearly by comparing 
eqs (4.2.11) and (4.4.10), which are repeated here (using t = nh) for emphasis. 

Analytical solution of ODE 


u =--■ c,(e A ‘Y7 * r 2 {e x - h ) n x 2 + ■■■ + c M {e XMh ) n x. M ^ [A] 7 (4.6.5a) 

Transient solution Steady state 

Numerical solution of ODE (one- root method) 

Un - c,(ai)"xi + c 2 {o 2 ) n X 2 + b c M (o M ) n x M + [A] 7 (4.6.5b) 

All of the c m ,x m ,/, and elements of j A ] in these eqs are identical. All of 
the one-root methods represented in Table 4.6.1 produce a result identical to that 
given by eq (4.6.5b) in which one inserts the appropriate o “signature” instead of eq 
(4.4.9). If a spurious root is generated, it adds another row to the expression for the 
transient with new constants, c m (which are again fixed by the initial conditions), 
but with the same eigenvectors and the same steady- state solution. 

4.7 Defective and/or Derogatory Matrices 

In general, the eigenvalues of a matrix may not be distinct, in which case the 
possibility exists that it cannot be diagonalized. If the eigenvalues of a matrix are 
not distinct, but all of the eigenvectors are linearly independent, the matrix is said 
to be derogatory, but it can still be diagonalized. A set of eigenvectors is linearly 
independent if a • x m + b ■ x n ^ x where m . ■£ n ^ k for any complex a and b 
and for all possible combinations of vectors in the set. However, if a matrix does 
not have a complete set of linearly independent eigenvectors, the matrix cannot be 
diagonalized, and it is said to be defective. A repeated root which causes a matrix 
to be defective will be referred to as a defective eigenvalue. Notice that, by this 
definition, a matrix can have an eigensystem with repeated eigenvalues that cause it 
to be both defective and derogatory. An example is given at the end of this section. 

Matrices which are not strictly diagonalizable can still be put into a compact 
form by a similarity transform, S, such that 

S~ 1 AS = \J] (4.7.1) 

and [J] is a Jordan matrix composed of blocks of submatrices spread along the 
diagonal, each submatrix having the form 
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A, 

0 


1 


0 



... 0 

At 1 : 

0 A, 0 

1 

••• 0 0 At 


(4.7.2) 


Use of this transformation is known as putting A into its Jordan canonical form. 
If a matrix A has two or more Jordan submatrices that have the same eigenvalue, 
the matrix is said to be derogatory. For each Jordan submatrix with an eigenvalue 
A i of multiplicity r, there exists one eigenvector. The other r — 1 vectors associated 
with this eigenvalue are referred to as principal vectors. 

In general, J has the form 


■j, [0] ••• [0] 

[o] J 2 : 

: 10] 

[0] ••• [0] J k 


(4.7.3) 


where there are at most k distinct eigenvalues. We use the term Jordan subblock, 
or simply Jordan block, to represent a matrix having the form given by eq (4.7.3) 
or to represent A, itself. For example, the matrix 


Ai 1 
Ai 1 

A 2 1 
A2 

^3 J 

[a. 


(4.7.4) 


27 



is both defective and derogatory, having: 

9 eigenvalues, 

4 distinct eigenvalues, 

6 Jordan blocks, 

6 linearly independent eigenvectors, 

2 principal vectors with Aj, 

1 principal vector with A 2 , 

2 defective eigenvalues, 

3 derogatory eigenvalues. 

Examples of defective eigenvalues occur in our subsequent development, see, for 
example, the discussion of the Gauss-Seidel method given later. 
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5. SOME PROPERTIES OF TRIDIAGONAL MATRICES 


5.1 Standard Eigensystem for Simple Tridiagonals 


In the following sections, tridiagonal banded matrices are prevalent. It is useful 
to list some of their properties. Many of these can be derived by solving the simple 
linear difference equations that arise in deriving recursion relations. 

Let us consider a simple tridiagonal matrix, i.e. a tridiagonal with constant 
scalar elements a, b, and c. If we examine the conditions under which the determi- 
nant of this matrix is zero, we find 


det[B(M:a,6,c)] = 0 if 

( 7Y17T \ 

— 1=0 , m = 1, 2, , M 

M + 1 / 


From this it follows at once that the eigenvalues of B(a,b,c) are 


A m = 6 + 2 Vac cos » 


m = 1,2 , • • • , M 


(5.1.1) 


(5.1.2) 


The right-hand eigenvector of B(a,b,c) that is associated with the eigenvalue A r 
satisfies the equation 


B(a, 6, c)x m \ m Xm 


(5.1.3) 


It is given by 


.7 - 1 

( a \~ . 

. / mn \ 

- sm 

j [ — ] 

\c / 

; \M + i ) _ 


m - 1,2 , • • • ,M (5.1.4) 


These vectors are the columns of the right-hand eigenvector matrix, the elements 
of which are given by 


X = 


{ a \ l ^ L • 

- sin 

jmn 

V c ) 

M + 1 


j = 1,2,--,M 
m — 1, 2, • • • , M 


(5.1.5) 


Notice that if a = —1 and c = 1, 
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The left-hand eigenvector matrix of B(a,b,c) can be written 


X 


- 1 


£) 


M + 

Notice that if a — — 1 and c = 1 


sin 


mj 7r 

M + 1 


m = 1,2, • • • , M 
j = 


(!) 


= e 


— t(m— 1 ) " 


5.2 Generalized Eigensystem for Simple Tridiagonals 

This system is defined as follows 


b c 
a b c 
a b 



*1 


« / 


x x 


*2 


d e f 


*2 


*3 

= A 

d e 


*3 




/ 




. _ 


d e 


Ki 


In this case one can show after some algebra that 

det[B(a - Ad, b - Ac, c - A/)j = 0 , if 
b X m e + 2\J{a A m d)(c A m /) cos ^ ^ ^ j ^ 0 ■> m 1? 2, ■ 


If we define 


e m = 


m 7r 

M + 1 


, Pm = COS 0 r 


(5.1.6) 

(5.1.7) 

(5.1.8) 

(5.2.1) 

! 

(5.2.2) 

(5.2.3) 
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eb - 2 {cd + a/)p^ ± 2 p m \[{ec - /6)(ea 6d) + ((ed - af)p, 


« 2 - 4 f d Pm 


(5.2.4) 


The right-hand eigenvectors are 


ci A ct 


C ~ Xmf 




sin [j'0 r 


m = 1,2, • • • ,M 

J = 1,2,--,M 


(5.2.5) 


These relations are useful in studying SOR methods later in this report. 

5.3 The Inverse of a Simple Tridiagonal 

The inverse of 5(a, 6, c) can also be written in analytic form. Let £>a* represent 
the determinant of 5(A/ : a, 6, c) 


= det(£?(M : a, 6, c)] (5.3.1) 

Defining Dq to be 1, it is simple to derive the first few determinants, thus 


D 0 = 1 
D i = 6 
Z ?2 = 6 2 — ac 
£>3 = 6 3 - 2a6c 

One can also find the recursion relation 


(5.3.2) 


Dm = bDsi-\ ~ o-cDm-i (5.3.3) 

Eq (5.3.3) is a linear OAE the solution of which was discussed in Section 4. Its 
characteristic polynomial P{E) is P(E 2 — bE + ac) and the two roots to P{o) = 0 
result in the solution 


Dm 


1 

f 

b + \/b 2 — 4ac 

Af+1 

6 — \/& 2 ~ 4ac 

M+l N 

> 

\Jb 2 — 4ac 

i 

2 


2 



M = 0, 1, 2, • • • 


(5.3.4) 
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where we have made use of the initial conditions Do = 1 and Dj = 6. In the limiting 
case when 6 2 - 4ac = 0, one can show that 


/ b \M 

Dm — + 1)( - j ; b 2 = 4 ac 

Then for M - 4 



Do 

— c£>2 

c 2 D\ 

- c z Do~ 

1 

- CfZ?2 

D iD? 

- cDi D\ 

c 2 D\ 

~Da 

a 2 D\ 

— aD\ D\ 

DiD\ 

- cD% 


- a 3 Do 

a 2 Di 

— aZ?2 

D 3 . 


and for M — 5 


D 4 

- cDz 

c 2 D 2 

- c 3 D\ 

c 4 D 0 

- aDz 

D j Dz 

~ cD\D2 

c 2 D x D 1 

- c 3 D, 

a 2 D 2 

aD j D2 

D2D2 

- CD2D] 

c 2 D 2 

- a 3 Di 

a 2 D\D\ 

- aD 2 D\ 

DzD\ 

- cDz 

a 4 Do 

- a 3 D\ 

a 2 D 2 

- &Dz 

d 4 


The general element d mn is 


Upper triangle: 

m=l,2, — 1 ; n = m+l,m + 2,--*,M 

dmn ~ D m -iDM-n{-c) n ”* / Dm 


Diagonal: 

n = m — 1,2,* • • , M 


dmm = Dm ^ 1 D m — m / Dm 


Lower triangle: 

m - n + l,n + 2,--,M ; n-l,2,--,M-l 

dmn — Dm - m D n _ ] ( — a) m n j Dm 


(5.3.5) 


(5.3.6) 


(5.3.7) 


(5.3.8) 
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5.4 Periodic or Circulant Tridiagonal Matrices 
Next consider the periodic tridiagonal shown in eq (1.1.4) 


B V (M : a,b,c) 


(5.4.1) 


The eigenvalues are 


X m 


b + (a 4- c) cos 




m = 0, 1 , 2, • • • , M - 1 


(5.4.2) 

Notice the slight shift in the index which makes the notation for the periodic analysis 
more convenient. The right-hand eigenvector that satisfies B p (a, b, c)x m = X m x m is 

= = , > = 0,1, 1 (5.4.3a) 

where i = >/—l. This can also be written 



j = - 1 (5.4.3b) 


The left-hand eigenvector matrix is 


X - 1 



m = 0, 1, • • • , M — 1 


(5.4.4) 


Notice the remarkable fact that the elements of the eigenvector matrices X and 
A' -1 for the periodic matrix do not depend on the elements a,b,c in the original 
matrix. In fact, all periodic (or circulant) matrices of order M have the same set 
of linearly independent eigenvectors. Further examination shows that the elements 
in these eigenvectors correspond to the elements in a complex harmonic analysis or 
complex discrete Fourier series. 

A full (or completely dense) circulant matrix of order M — 5 is shown in eq 
(5.4.5). 
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(5.4.5) 


B v 


bo 

b 4 

bz 

^2 

bi 


b i 62 63 64 

bn b \ b 2 63 

64 6n 6] ^2 

63 64 60 61 

62 ^3 ^4 60 


As we have just pointed out, the eigenvectors are given by eq (5.4.3) and are inde- 
pendent of the bj . However, the eigenvalues do depend on the matrix elements and 
they have the general form 


A 


m 


M l 

b J e l(2wjm/M) 

i=o 


(5.4.6) 


of which eq (5.4.2) is a special case. 


5.5 Special Cases Found From Symmetries 


Consider a mesh with an even number of interior points such as that shown in 
Fig 5.5.1a. One can seek from the tridiagonal matrix B(2M : a, 6, a, ) the eigenvector 
subset that has even symmetry. When spanning the interval 0 < x < n. For 
example, we seek the set of eigenvectors x m for which 


b 

a 






X\ 


T\ 

a 

b 

a 





*2 


I 2 


a 


* , 

a 



, 

— 

■ 




a 

b 

a 


*2 


X 2 





a 

b J 




_*1 . 


(5.5.1) 


This leads to the subsystem of order M which has the form 


B(M . ct , b , <r ) X Yn 


(5.5.2) 


where 
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b a 
a b 


a 


B(M : a, b, a) = 


• • a 

aba 
a b + a 


(5.5.3) 


By folding the known eigenvectors of B(2M : a, b, a) about the center, one can show 
from previous results that the eigenvalues of (5.5.3) are 


A 


m 


b + 2a cos 


/ (2m - 1 ) 7r \ 

\ 2M+1 ) 


m — 1 , 2 , • • • , M 


(5.5.4) 


and the corresponding eigenvectors are 


^ m 


= sin 


/ j (2m - 1 ) 7r \ 

V 2M 4 1 ) ’ 


J = 1 > 2, • • • , M 


(5.5.5) 


Imposing symmetry about the same interval 
but for a mesh with an odd number of points, 
see Fig. 5.5.1b, leads to the matrix 



Line of Symmetry 


0 0 0 

1 2 3 

1 2 3 

M 


o 

4 


o 


, 1 


5 


6 

M' 


B(M : a,b,a) 


b a 
aba 

a 

a 

aba 

2 a b 
(5.5.6) 


By folding the known eigenvalues of B(2M — 
1 : a, 6, a) about the center, one can show 
from previous results that the eigenvalues of 
(5.5.6) are 


a. An even-numb ered mesh. 

Line of Symmetry 

x = 0 

1 o o o o 

j’= 12 3 4 

j = 12 3 

M 

b. An odd-numbered mesh. 

Figure 5.5.1 - Symmetrical folds for 
special cases. 


. 1 

5 

M' 
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Am = b + 2a cos 


/ (2m — 1 ) 7r \ 

V 2 M / 


, m = 1,2, • • • ,M 


(5.5.7) 


and the corresponding eigenvectors are 


Xm - sin 


.... ( li 2m ~ 1)* \ ■ _ ! O 

v 2M j ’ 




(5.5.8) 


5.6 Special Cases Involving Boundary Conditions 

We consider three special cases for the matrix operator representing the 3-point 
central difference approximation for the second derivative d 2 /dx 2 at all points away 
from the boundaries, combined with special conditions imposed at the boundaries. 


Note: In every case 

m 

J 

— 2 + 2 cos (a) 


= 1,2, 

= 1,2 ,--,M 

= -4 sin 2 (a/2) 


W'hen the boundary conditions are Dirichlet on both sides, 


- 2 1 
1-2 1 
1-2 1 
1-2 1 
1 - 2 



(5.6.1) 


When one boundary condition is Dirichlet and one is (first-order) Neumann 


2 1 
1-2 1 
1-2 1 
1-2 1 
1 - 1 


(2m - 1) 

2M~ 


jc jyi — sin 


/ (2m - 1 ) 7T \ 

V 2 aT+T j 


\ 


> 


(5.6.2) 
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When the boundary conditions are Neumann on both sides 


1 1 
1-2 1 
1-2 1 
1-2 1 
1 - 1 


-2 + 2 cos 


(m - 1 )?r 
M 


Xm = COS 




(m - l)7r\ 

M ) 


Notice that only the last matrix is singular. 


37 



6. CLASSICAL RELAXATION 


6.1 The Converged Solution, the Residual, and the Error 

Most standard texts, see refs. 4,5,6, approach the subject of relaxation by 
considering what we have defined as the preconditioned form, see eq (2.2.1) and 
(2.3.2), 


A4>- f = 0 


( 6 . 1 . 1 ) 


and applying to it the iterative process 


t n+1 = \I+MA}i n -Mf (6.1.2) 

In the terminology of Section 3, M is the secondary conditioning matrix equal to 
H -1 in eq (3.3.1). Equation (6.1.2) is usually rewritten in the form 


k+l = GK - #71 

where G = / + HA J 


(6.1.3) 


In this notation, G is the basic iteration matrix and its eigenvalues, which we des- 
ignate as <7 m , determine the convergence rate of a method. The converged solution 
to all three eqs (6.1.1), (6.1.2), and (6.1.3) is 


ico = A~'f (6-1-4) 

The error at the nth iteration is defined as 


= A~'f (6.1.5) 

The residual at the nth iteration is defined as 

r n = A4> n -f (6.1.6) 

Multiply eq (6.1.5) by A from the left, and use the definition in eq (6.1.6). There 
results the relation between the error and the residual 


Ae n - r„ - 0 


(6.1.7) 
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It is not difficult to show that 


e n+ i^(7e„ (6.1.8) 

In all of the above, we have considered only what are usually referred to as 
stationary processes. We are in fact much more interested in nonstationary pro- 
cesses. but our approach to them is not standard so we defer the discussion of these 
to Section 7. 


6.2 Point Operator Schemes in One Dimension 

Let us consider three classical relaxation procedures for the one dimensional 
equation 


0= 9 (x) (e.2,) 

which, with three-point central differencing and Dirichlet boundary conditions, re- 
duces to the model equation 


B(l,-2,1)4>= f (6.2.2) 

where / ~ A x 2 g. These methods are very well known but the terminology is not 
universal. For example, the Gauss-Seidel method is sometimes called the Lieb- 
rnan method and the Point-Jacobi method has been referred to as the Richardson 
method. 

What we refer to as the Point-Jacobi method is expressed in point operator 
form for the one dimensional case as 


u 


n -b 1 ) 


(n) , In) _ , 

y-i + j+i 


(6.2.3) 


The Gauss-Seidel method is 


u 


(n4 1) 


,(”+!) , ,» 

S-i 


+ ~ u 


(6.2.4) 


The method of successive overrelaxation (SOR) is usually expressed in two steps as 
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(6.2.5a) 


u, = - 
3 2 


u 


(n-j 1 ) (n) 

u] - w] -I w 


(n+1) (n) _ , 

j - 1 + u y + i h 

(n) 


U,' U 
J j 


but it can also be written in the single line 


2 ”j-* ■ ' "J 2 J ' +1 2 

6.3 The Convergence Rates 


(6.2.5b) 


The usual measure of the convergence rate is the eigenvalue cr m of G, see eq 
(6.1.3), having maximum absolute value. Thus 


Convergence 


\o r 


; m = 1, 2, • • • , M 


(6.3.1) 


These values are well known for Laplace’s equation using the three methods just 
defined. They are 


Point- Jacobi 


Gauss-Seidel 


SOR 


7T 

= COS . . , 

m <max ' M + 1 


\fT ! — 

( * Y 

m Imax 

Cuo 1 _ ’ 1 

\M + 1 / 


l 2 


m I max 


— w opt J 


W 0f .t = 2/ 


1 + sin 


7T 


AJ + 1 


(6.3.2a) 


(6.3.2b) 


(6.3.2c) 


41 



7. ODE APPROACH TO CLASSICAL RELAXATION 


7.1 ODE Form of the Classical Methods 

The three iterative procedures defined by eqs (6.2.3), (6.2.4), and (6.2.5) obey 
no apparent pattern except that they are easy to implement in a computer code 
since all of the data required to update the value of one point are explicitly available 
at the time of the update. Now let us study these methods as subsets of ODE as 
formulated in Section 3. Insert the model equation (6.2.2) into the ODE form 
(3.3.2). Then 


H— = B(l, -2, \)<f> — f 


(7.1.1) 


As a start, let us use for the numerical integration the explicit Euler method 


<^n+ 1 — <t>n + h(f)' n 


(7.1.2) 


with a step size, /i, equal to 1. We arrive at 


H(k+,-K) = B(l,-2,l)i n -f (7.1.3) 

It is clear that the best choice of H from the point of view of matrix algebra is 
B( 1,— 2,1) since then multiplication from the left by — B -1 (l,— 2, 1) gives the 
correct answer in one step. However, this is not in the spirit of our study, since 
multiplication by the inverse amounts to solving the problem by direct methods 
without iteration. The constraint on H that is in keeping with the formulation of 
the three methods described in Section 6 is that all the elements above the diagonal 
(or below the diagonal if the sweeps are from right to left) are zero. If we impose 
this constraint and further restrict ourselves to banded tridiagonals with a single 
constant in each band, we are led to 

0)(^ +1 - 4> n ) = B(l, —2, l)(j> n - f (7.1.4) 

where /? and w are arbitrary. With this choice of notation the three methods 
presented in Section 6 can be identified using the entries in Table 7.1.1. 
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TABLE 7.1.1 - VALUES OF (3 and u IN EQ (7.1.4) THAT 
LEAD TO CLASSICAL RELAXATION METHODS 

/ 3 u: 

Method 

Equation 

0 1 

Point- Jacobi 

6.2.3 

1 1 

Gauss- Seidel 

6.2.4 

1 2/[l + sm ( M+1 )] 

Optimum SOR 

6.2.5 


The fact that the values in the tables lead to the methods indicated can be 
verified by simple algebraic manipulation. However, our purpose is to examine the 
whole procedure as a special subset of the theory of ordinary differential equations. 
In this light, the three methods are all contained in the set of ODE 

~ m B- 1 (-y3,^,O)[B(l,-2,l)0-/] (7.1.5) 

and appear from it in the special case when the explicit Euler method is used for 
its numerical integration. The point operator that results from the use of the Euler 
scheme, 4> n -\-\ — + h(p' n , is 




(n+l) _ 


o 


+ 


(jj 

-( h 
2 ' 


0)^-i ) - (M- i)^- n) ) + (y^+i) 


u )h 

T ’ 1 

(7.1.6) 


This represents a generalization of classical relaxation that results from the numer- 
ical solution of ODE. 


7.2 The A Eigenvalues and the Error 

It is at this point that we start to deviate from the usual presentation of relax- 
ation. In our approach it is essentia] that we first identify the A eigenvalues of the 
iterative procedure, and then determine the o eigenvalues as a function of A h. The 
A eigenvalues are fixed by the basic matrix in eq (2.1.1), the preconditioning matrix 
in (2.2.1), and the secondary conditioning matrix in (3.3.2). The o eigenvalues are 
fixed for a given A h by the choice of numerical method as illustrated in Table 4.6.1. 
Throughout the remaining discussion we will refer to the independent variable t as 
“time”, even though no true time accuracy is necessarily involved. 


44 



The general form of the relaxation procedure, as we are reviewing it, is 


^ = H~' 
dt 


C 


A b 4> — f i, 


H~ l \A4> - f) 


(7.2.1) 


From Section 4, we see that, if the eigenvectors of \H 1 A] are linearly independent, 
the solution can be written as 


4> = cie Xlt Xi + h CMe Awt XM +0oo (7.2.2) 

s. — V~ — ^ 

error 

where what is referred to in time-accurate analysis as the transient solution (see 
eq (4.2.11)), is now referred to in relaxation analysis as the error (see Section 6). 
We see that the problem in the relaxation subset of ODE is to remove the transient 
solution from the general solution in the most efficient way possible. 

7.3 Stationary Processes 

If H and C in eq (7.2.1) are independent of t, that is, are not changed through- 
out the iteration process, the method is referred to, in relaxation terminology, as 
stationary. The generalization of this in our approach is to make h, the “time” step, 
a constant for the entire iteration. 

Suppose the explicit Euler method is used for the time integration. Then, from 
Section 6 and eq (7.2.2), the numerical solution can be expressed, after n steps, as 


4> n — c lXi(l + A i /z) + • • • + c m x m (l + X m h) + • • • + CJW^m(1 + X^fh) +</>oo 

' V 

error 

(7.3.1) 

The initial amplitudes of the eigenvectors are given by the magnitudes of the 
c m . These were fixed by the initial guess. In general it is assumed that any or all 
of the eigenvectors could have been given an equally “bad” excitation by the initial 
guess, so that we must devise a way to remove them all from the general solution 
on an equal basis. Assuming that [H~ 1 A] has been chosen (that is, an iteration 
process has been decided upon), the only free choice remaining to accelerate the 
removal of the error terms is the choice of h. The three methods represented in 
Table 7.1.1 have all been conditioned by the choice of H to have an optimum h 
equal to 1 for a stationary iteration process. 
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7.4 Nonstationary Processes 


In classical terminology a method is said to be nonstationary if the conditioning 
matrices, | H 'C], are varied at each time step. This would not change the steady 

state solution but it can greatly affect the convergence rate. In our ODE 

approach this could also be considered and would lead to a study of equations with 
nonconstant coefficients. It is much simpler, however, to study the case of fixed 
\H " 1 C ] but variable step size, h. This process changes the Point-Jacobi method 
to Richardson’s method in standard terminology, see ref. 4. For the Gauss-Seidel 
and SOR methods it leads to processes that are superior to the stationary methods 
but, to the author’s knowledge, unpublished. 

The nonstationary form of eq (7.3.1) is 

N N 

<j> n = C\X\ 1 1 (l -f Aj/t n ) + • • • + c m i m J | (1 + X m h n ) 

n --- 1 n~ 1 

N 

+ • • ’ + C M X M n (1 + Aj \ih n ) + <i> 

oo 

n— 1 


(7.4.1) 


where the symbol n stands for product. Since h n can now be changed at each step, 
the error term can theoretically be completely eliminated in M steps by taking 
= -1/A m , for m = 1,2 This concept of eigenvector annihilation is 
discussed in Section 8. 

Let us consider the very important case when all of the X m are rea l and negative 
(remember that they arise from a conditioned matrix so this constraint is not un- 
realistic for quite practical cases). Consider one of the error terms taken from 

M N 

= 4*n ^oo = ^ ^ x m || (l “I” A m h n ) (7.4.2) 

m— 1 n= 1 

and write it in the form 

N 

Cm%m P e (^m) — Jj| (l d - ^m^n) (7.4.3) 

n- 1 

where P € signifies an u Euler” polynomial. Now focus attention on the polynomial 
(Pe) N W = (1 + ^ 1 A)(1 +h 2 \)---{l + h N \) (7.4.4) 
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treating it as a continuous function of the independent variable A. In the annihi- 
lation process mentioned after eq (7.4.1) we considered making the error exactly 
zero by taking advantage of some knowledge about the discrete values of A m for a 
particular case. Now we pose a less demanding problem. Let us choose the h n so 
that the maximum value of (P € )n(\) is as small as possible for all A lying between 
A n and Aj, such that Af, < A < A a < 0. Mathematically stated, we seek 


max |(P e )^r(A)| = minimum ; with(P e )^ (0) = 1 

A /, < A < A (l 

This problem has a well known solution due to Markov. It is 


{ P e)»{ A) 




(7.4.5) 


(7.4.6) 


where 


T N (y) = cos(N arccosy) 

are the Chebyshev polynomials along the interval — 1 < y < 1 and 


(7.4.7a) 


r w (w) = | ' 


■ — 1 N 


y + vV - ij +^[y-\/y^ 


-1 N 


(7.4.7b) 


are the Chebyshev polynomials for |y| > 1. In relaxation terminology this is gener- 
ally referred to as Richardson’s method, see ref. 6, and it leads to the nonstationary 
step size choice given by 


= ^(-A t -A, + (A t -A,)eos }■ »=1,2,--AT (7.4.8) 


Remember that all A are negative real numbers representing the magnitudes of A m 
in an eigenvalue spectrum. 

The error in the relaxation process represented by eq (7.4.1) is expressed in 
terms of a set of eigenvectors, x m , amplified by the coefficients c m f](l + A m h n ). 
With each eigenvector there is a corresponding eigenvalue. Equation (7.4.8) gives 
us the best choice of a series of h n that will minimize the amplitude of the error 
carried in the eigenvectors associated with the eigenvalues between A^ and A n . 
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As an example for the use of eq (7.4.8), let us consider the following problem: 


Minimize the maximum error asso- 
ciated with the A eigenvalues in the 
interval -2 < A < — 1 using only 3 
iterations. 


(7.4.9) 


The three values of h which satisfy this prob- 
lem are 



and the amplitude of the eigenvector is re- 
duced to 

(P e ) 3 (A) = r 3 ( 2 A + 3)/r 3 (3) (7.4.11) 

where 

7s ( 3) = {>3 + \ r 8\ 3 + [3 - v / 8J 3 }/2 « 99 

(7.4.12) 

A plot of eq (7.4.12) is given in Fig. 7.4.1 and 
we see that the amplitudes of all the eigen- 
vectors associated with the eigenvalues in the 
range -2 < A < -1 have been reduced to less 
than about 1% of their initial values. 

Return now to eq (7.4.1). This was derived 
from eq (7.2.2) on the condition that the ex- 
plicit Euler method, eq (7.1.2), was used to 
integrate the basic ODE. If instead the im- 
plicit trapezoidal rule 


h x = 4/(6- V3) 
h 2 = 4/(6 - 0) 
/13 = 4/(6 + v^3) 



Figure 7.4.1 - Richardson’s 
method for 3 steps, 
minimized for -2 < A < — 1. 


<£n+l — 4>n + + 


(7.4.13) 


is used, the nonstationary formula 


<t>N 


M 

771=1 


N 

17 1 

f 1 — 

2 bn ^ 

11 1 

n - 1 

U- 

2 ) 


+ ^OO 


(7.4.14) 
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would result. This calls for a study of the rational “trapezoidal” polynomial, Pt, 


(Pi)n( A) = 


n P + 

“A 1 l h ^) 


(7.4.15) 


under the same constraints as before, namely that 


max |(.Pi)jv(A)| = minimum ; 

A i, ^ A A „ 

with (P/)/v(0) - 1 


(7.4.16) 


The optimum values of h can also be found for 
this problem, see ref. 6, but we settle here for 
the approximation suggested by Wachspress 




(»-l)/(JV-l) 


n= 1,2, •,7V 

(7.4.17) 


This process is also applied to problem (7.4.9). 
The results for (P<) 3 (A) are shown in Fig. 7.4.2. 
The error amplitude is about 1 /5 of that found 
for (P e ) 3 (A) in the same interval of A. 


h i — 1 

h 2 = V 2 
h z = 2 



Figure 7.4.2 - Wachpress 
method for 3 steps, 
minimized for -2 < A < — 1. 


7.5 Eigensystems of the Classical Methods 

Before we carry the ODE analysis further, it is instructive to inspect the eigen- 
vectors and eigenvalues in the \H~ X A\ matrix for the three classical methods repre- 
sented by eqs (6.2.3), (6.2.4), and (6.2.5). Using eq (7.1.5) we see that this amounts 
to solving the general eigenvalue problem 

P A m ^4x w (7.5.1) 
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for the special case 


B(-f3, - ,0)x m = -2, l)x m (7.5.2) 

id 

Eq (7.5.2) is a special case of eq (5.2.1) so the solution is given in Section 5. The 
three special cases are considered below. To illustrate the behavior we take M = 5 
for the matrix order. This special case makes the general result quite clear. 

7.6 The Point-Jacobi System 

If (3 = 0 and w = 1 in eq (7.1.5), the ODE matrix reduces to simply B(|, —1, |). 
The eigensystem for this matrix is given in eqs (5.1.1) through (5.1.8). One finds 
the following. 


Columns of eigenvectors 


1/2 

>/3/2 

1 

v/3/2 

1/2 

v3/2 

V' 3/2 

0 

- v^/2 

- v/3/2 

1 

0 

- 1 

0 

1 

V3/2 

v / 3/2 

0 

V / 3/2 

- y / 3/2 

1/2 

- v / '3/2 

1 

v/3/2 

1/2 


Eigenvector number 


1 


2 3 4 


5 


(7.6.1a) 


Eigenvalue 


-2 + \/ 3 ] -1 


§[-2+v/3| 


The left hand eigenvector set that is consistent with this normalization is 


Eigenvector number 


1/2 

v/3/2 

1 

v/ 3/2 

1/2" 

} 

1 

V3/2 

v / 3/2 

0 

- V3/2 

- v/3/2 

} 

2 

1 

0 

- 1 

0 

1 

} 

3 

v/3/2 

- v/3/2 

0 

V3/2 

- v/3/2 

} 

4 

. 1/2 

- v/3/2 

1 

- v/3/2 

l/ 2 . 

} 

5 


(7.6.1b) 
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If we inspect eq (7.4.1), we notice that the initial error content is given by 


*0 = <t > 0 - oo = X( 

Insert eq (5.1.5) for X, then 


(7.6.2a) 


M 

(eo )j = c msin 

m— 1 



j = 1»2,- -,M 


(7.6.2b) 


Tliis simply states that the elements in 
the initial error vector, or the initial tran- 
sient terms in the ODE, are given by a 
sine transform of the amplitudes of the 
eigenvector content. Similarly 



represents c as a sine synthesis of e a . This 
is a very “well-behaved” eigensystem with 
linearly independent eigenvectors and dis- 
tinct eigenvalues. The first 5 eigenvectors, 
which are simple sine waves, are shown in 
Fig. 7.6.1a. The eigenvalues are given by the 
equation 



, / mn \ 

A m = -l+cos^— — j (7.6.4) 

The functional dependence of a on A for the explicit Euler method is < 7 m = 1 + A m h. 
Thus the a, A relation can be plotted for any h. The plot for h = 1, the optimum 
stationary case, is shown in Fig. 7.6.1b. 
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7.7 The Gauss-Seidel System 


If /? and u> are equal to 1 in eq (7.1.5), the matrix eigensystem evolves from 
the relation 

B(l, -2, l)£ m = -2,0 )x m (7.7.1) 

which can be studied using the results in Section 5.2. One can show that the 
\HA~ l ] matrix for the Gauss-Seidel method, , is 

[A\gs - B-'(-l,2,0)B(l,-2,l) = 

1 
2 

3 (7.7.2) 

4 

5 

M 

The eigenvector structure of the Gauss-Seidel ODE matrix is quite interesting. If M 
is odd there are (At + l)/2 distinct eigenvalues with corresponding linearly indepen- 
dent eigenvectors, and there are (M - l)/2 defective eigenvalues with corresponding 
principal vectors. The Jordan canonical form for M =■ 5 is shown in Fig. 7.7.1. 




Figure 7.7.1 - Jordan canonical form 



- 1 1/2 

0 -3/4 1/2 

0 1/8 -3/4 1/2 

0 1/16 1/8 -3/4 1/2 

0 1/32 1/16 1/8 -3/4 1/2 

0 1 / 2 m 
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The right-hand column vector set, X, is shown in eq (7.7.3), where the * stands 
for a principal vector in the defective system, and the corresponding left hand row 
vectors, X -1 , are shown in eq (7.7.4). Neither set has been normalized in any 
rational way. 



16 


16 1 

0 

0 



24 


8 0 

2 

0 


X = 

24 


0 0 

- 1 

4 



18 


- 2 0 

0 

-4 



9 


- 1 0 

0 

1 


Vector number 

1 

✓ 

2 3 


* 

* 

Eigenvalue 

- 1/4 

-3/4 - 

1 

- 1 

- 1 



'o 

1/144 

1/72 

1/54 

1/54 



0 

1/16 

1/8 

0 

-1/2 

x- 1 

— 

1 

- 10/9 - 

20/9 

- 8/27 

208/27 



0 

5/30 

- 2/3 

- 2/9 

16/9 



0 

0 

0 

-1/6 

1/3 


(7.7.3) 


(7.7.4) 


The equation for the nondefective eigenvalues in the ODE matrix is (for odd M) 

, o, m-n . M + 1 . , 

A m = -1 +cos 2 (^-j) ; m = 1, 2, • • • , — - — (7.7.5) 

and the corresponding eigenvectors are given by 


X 


m 


COS 


m 7r 

M + 1 


• 3 - 1 

sin 



m — 1 , 2 , 


M+ 1 


(7.7.6) 


The eigenvectors are quite unlike the Point-Jacobi set. They are no longer 
symmetrical. They produce waves that are higher in amplitude on one side (the 
updated side) than they are on the other. Fig. 7.7.2a shows the eigenvectors corre- 
sponding to the lowest three eigenvalues for M = 11, 23, and 47. Notice that they 
do not represent a common family for different values of M . The o m produced from 
the A m by the explicit Euler method varies with h. The a, A relationship for h. — 1, 
the optimum stationary case, is shown in Fig. 7.7.2b. 
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7.8 The SOR System 


If j) - 1 and 2 ,/u? - x in eq (7.1.5), the ODE matrix is ti " 1 (- 1, x,0)R(l, -2, 1). 
One can show that this can be written in the form given below for M — 5. The 
generalization to any M is fairly clear. The H 1 A matrix for the SOR method, 
A\ soh = R-'( -l,x, 0)5(1, -2,1), is 


- 2x 4 x 4 0 0 0 

- 2x 3 + x 4 x 3 — 2x 4 x 4 0 0 

- 2x 2 + x 3 x 2 - 2x 3 + x 4 x 3 - 2x 4 x 4 0 

- 2x + x 2 x - 2x 2 + x 3 x 2 - 2x 3 -f x 4 x 3 - 2x 4 x 4 

-2 + i 1 - 2x + x 2 x - 2x 2 + x 3 x 2 - 2x 3 + x 4 x 3 - 2x 4 

(7.8.1) 


Eigenvalues of the system are given by 


Am - ~ 1 + 


Upm + Zm 
2 


z m —[4(1 — <*?) + W 2 p^,j 1/2 

p m = cos[m7r /(M + 1)] 



(7.8.2) 
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If u = 1, the system is Gauss-Seidel. If 4(1 — w) +u> 2 p m <0 ,z m and A m are complex. 
If u is chosen such that 4(1 - u) + u; 2 p 2 = 0,w is optimum for the stationary case, 
and the following conditions hold 

1) Two eigenvalues are real, equal and defective 

2) If M is even, the remaining eigenvalues are complex and occur in conjugate 
pairs 

3) If M is odd, one of the remaining eigenvalues is real and the others are complex 
occurring in conjugate pairs. 

One can easily show that the optimum w for the stationary case is 


u op t = 2/ 


1 + sin 



(7.8.3) 


and for u> = w G pt 


A 


m 



- 1 


Xm =Cm 1 sin 



where 


Cm — 


k , jypt 


Pm + i\f?\ 


Pm 

(7.8 


4) 


If the explicit Euler method is used 
to integrate the ODE, o m = 1 - 
h + h and if h -- 1, the optimum 
value for the stationary case, the a, A 
relation reduces to that shown in 
Fig. 7.8.1. This illustrates the well 
known fact that for optimum sta- 
tionary SOR all the \o m \ are iden- 
tical and equal to u opi - 1. 



Figure 7.8.1 - The a, A relation for 
optimum stationary SOR, M — 5, h = 1. 
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7.9 Solution of the ODE Forms of the Classical Methods 


The ODE approach to relaxation can be summarized as follows. The basic 
equation to be solved came from some time accurate derivation 


A b u-f h = 0 (7.9.1) 

This equation is preconditioned in some manner which has the effect of multiplica- 
tion by a conditioning matrix C giving 


A4>-f- 0 


(7.9.2) 


An iterative scheme is developed to finding the converged, or steady-state, solution 
of the set of ODE 


H 


d(j> 

dt 


= A 4> - f 


(7.9.3) 


This solution has the analytic form 


Transient term or I f Steady state term, 1 (7 9 4) 

error, e r ) \ ^ = A -1 / J 

The three classical methods, Point-Jacobi, Gauss-Seidel, and SOR, are identified 
for the one-dimensional case by eq (7.1.5) and Table 7.1.1. Their solution for M = 5 
is written below. For all cases involving linear systems, higher values of M, or higher 
dimensions, do not alter the fundamental nature of the solutions. 

P oint-Jacobi : The explicit Euler method is used to solve 



d4> 

dt 



2,1 )j-f 


(7.9.5) 


The eigenvectors are all real and linearly independent, see eq (7.6.1). For M — 5 
they can be written 
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' 1/2 ' 


y/3/2 


r 


a/3/2 ' 

v/3/2 


v/3/2 


0 


- \/3/2 

1 

,x 2 = 

0 

1 X 3 = 

- 1 

,x 4 — 

0 

\/3/2 


- v/3/2 


0 


v/3/2 

. ! / 2 . 


. -\/3/ 2. 


1 


_ - v/3/2. 


X 5 


1/2 

- V3/2 

1 

- V3/2 

1/2 


The corresponding eigenvalues are, from eq (7.6.4) 

A, = -1 + & = —0.134 • • • ' 
A 2 = 1 + 5 = 0.5 

A 3 = - 1 =-1.0 > 

A 4 =-l-± =-1.5 

A 5 = - 1 - ^ = -1.866- •• , 

The numerical solution written in full is 


t n - too = cj|l - (1 - —)h] n x | 

+ call - (1 - IW'ii 

+ C 3 !l - (1 )fc] n X 3 
+ c 4 [l - (1 + -)h} n x 4 

+ ^5 i 1 ~ (1 + “J^pXs 

Gauss-Seidel : The explicit Euler method is used to solve 

^ = B- I (-l,2,0)[B(l,-2,l)^-7] 


(7.9.6) 


(7.9.7) 


(7.9.8) 


(7.9.9) 
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The eigenvectors and principal vectors are all real. For odd M, (M + l)/2 of them 
are linearly independent eigenvectors and (M + l)/2 are principal vectors. For 
M 5 they can be written, see eq (7.7.6) 


1/2 


x/3/2 


V 


0 " 


0 " 

3/4 


V^/4 


0 


2 


0 

3/4 

,12 = 

0 

,*3 = 

0 

, X 4 -- 

- 1 

,*5 = 

4 

9/16 


- a/3/16 


0 


0 


- 4 

9/32 


- \/3/32_ 


0 


0 


1 


(7.9.10) 

The corresponding eigenvalues are, from eq (7.7.5) 

A ! - - - 1/4 = -0.25 
A 2 = - 3/4 = -0.75 

A 3 = — 1 =-1.00 (7.9.11) 

(4) \ Defective, linked to 

(5) J A 3 Jordan block 

The numerical solution written in full is 


4>n - <^oo 

. 3 /i. - 

■+ C 2 (l - — ) x 2 


cz( 1 -h) n + c 4 h^{ 1 -fe)"- 1 +c 6 fc 


n(n - 1) 
~ 2! 


(i - />) 


n-2 


X3 


-j- c 4 (l - h) n T C5/1 — (1 — ll) 
+ 05(1 - h) n X 5 


to — 1 


X4 


(7.9.12) 


See Section 8.2 for an interesting experiment with this defective system. 
SOR : The explicit Euler method is used to solve 


dd> 

dt 



B(l, - 2, !)</> - / 


where cjq — 2/ { 1 + sin[7r/(M + ])]} 


(7.9.13) 
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For odd M, three of the vectors are real, two of these are eigenvectors and one 
is a principal vector. The remaining vectors are all complex, linearly independent 
eigenvectors. For M — 5 they can be written, see eq (7.8.4) 


1/2 


’ - 6 


>/3(l )/2 


1 

1/2 


9 


v/3(l± *V2)/6 


0 

1/3 

9 *^2 

16 

> x 3,4 — 

0 

9 

1/3 

1/6 


13 


V3(5± t'\/2)/54 


0 

1/18 


6 


>/3(7 ±4iy/2)/ 162 


.1/9. 


(7.9.14) 

The corresponding eigenvalues are 


A, =r -2/3 

( 2 ) Defective linked to Aj 

A 3 =-(10-2 v / 2 j)/9 (7.9.15) 

A 4 = (10 l 2V2*)/9 

A 5 = -4/3 

The numerical solution written in full is 


4>n-4 > oo = [Cl (1 - 2h/Z) n + c 2 nh( 1 - 2h/3) n ~ 1 )x i 
+ c 2 {l-2h/3) n i 2 

+ c 3 [l - (10 - 2\/2i)h/9] n X3 (7.9.16) 

+ c 4 [l - (10 -t- 2v / 2*)^/9] ri ^ 4 
+ 05(1 - 4h/3) n x$ 
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8. EIGENVECTOR ANNIHILATION 


8.1 Introduction 

Let us consider a conditioned matrix that has a set of linearly independent 
eigenvectors. We have seen that the ODE approach to relaxation first casts the 
process into a form which has the exact analytical solution 


4> n = c 1 (e x ' h ) n i 1 + ... + Cm (e A ’" h )"i m + ... + c M (e AA,fc ) n ^M + <L (8.1.1) 

and then awaits the choice of a numerical integration method to relate the A eigen- 
values of the differential equation to the o eigenvalues of the finite difference equa- 
tions. If there is one a for each A, we find 

M ( N \ 

fin ~ ^ ^ { c m}s 1 1 cr m(^m^n) \ x m + ^oo (8.1.2) 

m— 1 l n - 1 J 

where we have taken advantage of the fact that the step size h can be varied from 
step to step. Clearly the problem is to remove the eigenvector content of the tran- 
sient solution which is identical to the error of the relaxation process. Three ways 
of doing this are immediately evident: 


(1) Make a good initial guess, that is, make c m « 0. 

(2) Make all jcrj < 1 and iterate many times; that is, make N large. 

(3) Try to choose h n so that |<r m | ss 0 aL least once in the process. 

Of these, (1) is obvious if it is possible, (2) is the typical approach of stationary 
methods, and (3) is the basis for nonstationary methods. 

In addition to the above, there are two more subtle ways for diminishing the 
error: 

(4) Mixing the eigenstructure by applying appropriate conditioning ma- 
trices and restarting the process. This makes a sophisticated use of (1) 
above. 

(5) Using multiple grids to surface and destroy the eigenvectors associated 
with the lower space-frequencies. 
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These two concepts are discussed in Sections 9 and 10. 

8.2 Selective Eigenvector Annihilation 

One of the most direct ways to eliminate the eigenvector content of the error 
in a relaxation scheme is to set h = -1/A m for every A m eigenvalue in the system 
and iterate M steps. This is an extremely poor strategy in practical application 
but it is worth discussing from a theoretical viewpoint. It is of course, in principal, 
a direct solution. 

Consider the Point-Jacobi solution given by eq (7.9.8). In the 5-point mesh, 
5 sweeps through the mesh with h equal to 7.464. . ., 2, 1, 0.666. . 0.536. . . would 
annihilate all 5 of the eigenvectors. In carrying out the calculation, the worst situ- 
ation from computer considerations would occur for the fifth eigenvector. Its initial 
amplitude would be multiplied consecutively by ( 12.93 . . .) ( -2.73 . . .) ( .866 . . .) 
(-.245 . . .) (0). In this trivial case the first tw'o multiplies would not cause trouble, 
but in application to large systems such a strategy would be highly unstable because 
computer hardware limits the accuracy to which numbers can be represented. No- 
tice. however, that the eigenvectors associated with the highest (M - l)/2 eigenval- 
ues can all be annihilated in this manner without amplifying any_vector at any step . 
Only when annihilating the eigenvectors associated with the lowest (M - l)/2 eigen- 
values are the amplitudes of other vectors amplified. 

Consider next the Gauss-Seidel solution given by eq (7.9.12). In an M point 
mesh, M sweeps are still required to annihilate all the possible error vectors. If M 
is odd, (M + l) / 2 of the sweeps with h — 1 are required to remove the principal vec- 
tors associated with the defective eigenvalues. The remaining annihilations require 
(M - l)/2 sweeps with h — -1/A m . One can observe the initial growth of the 
principal vectors for h i / 1 caused by the factor h k (n)(n — 1) • • • (n - k + 1 )/fc! 
by performing simple numerical experiments. Since the eigenvalues of the Gauss - 
Seidel system are all real and there are only half as many of them as there are in 
Point-Jacobi, one might question the second sentence in this paragraph and be led 
to the conclusion that the reduction in the number of eigenvalues could be used 
to advantage in annihilation. That such is not the case can be demonstrated by 
using the simple program written in BASIC and presented in Fig. 8.2.1. Recall the 
system of eigenvectors and principal vectors give by eq (7.9.10) and consider the 
pathological nature of eq (7.9.12) when h = 1 and the exponent of (1 - h) is zero. 
Run the BASIC program using — j0, 0,4, 4. lj as input. Notice that after one 

iteration the second principal vector x 4 appears. A second iteration produces the 
eigenvector x 3 and a third j iteration is required to annihilate the entire 

Jordan block. 
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10 DIM U(6) 

20 PRINT “ENTER INITIAL VALUES” 

30 INPUT U(1),U(2),U(3),U(4),U(5) 

40 FOR J - 1 TO 5 

50 IJ(J) = ,5*(U(J-1) + U(J + 1)) 

60 NEXT J 

70 PRINT U(1),U(2),U(3),U(4),U(5) 

80 PRINT 

90 INPUT “TYPE Y TO ITERATE AGAIN, OTHERWISE NEW INI- 
TIAL V” ;CH$ 

100 IF CHS = “Y” GOTO 40 
110 GOTO 20 
120 END 


Figure 8.2.1 - BASIC program to illustrate defective nature of Gauss-Seidel 
eigensystem. 

Input 0,0,4 ,-4,1 -► 0,2,-1,0,0-^1,0,0,0,0-^Q 


Finally, consider the SOR solution given by eq (7.9.16). It is again true that 
M sweeps will annihilate all of the possible error vectors, only now one must use 
complex values for the step size h. In present day computers this requires more 
arithmetic and storage. Nevertheless, it may have features that make it attractive 
for practical use. This is being investigated. 

We now make some observations. First of all, it is usually unnecessary to re- 
move completely the error from the general solution. We consider it satisfactory to 
reduce the error below some threshold. Second, the exact values of A m are generally 
unknown and not practical to compute. It is therefore unnecessary and impractical 
to set h = -1/A m for m = 1,2, , A/; a few well chosen h ' s will greatly reduce the 
amplitudes of whole clusters of eigenvectors associated with nearby A’s in the A m 
spectrum. For example, in the Point-Jacobi case on the model problem using only 
3 h ' s that are given by the — 1 / A m at the two ends and in the middle of the highest 
half of the |A m | range, reduces the amplitudes of all the eigenvectors associated with 
that range to about 1% (or less) of their initial values. The amplitudes of all the 
other eigenvectors coupled into the system are also reduced, but some by very little. 
Fortunately, this phenomenon is not sensitive to the particular choice of h. The val- 
ues h — ^ , | , 1 are about as good as the Chebyshev choice h = 0.517 . . . , |, 0.937 . . . 
given in Fig. 7.4.1. These observations are much more significant when extended 
to 2 and 3 dimensions. In these cases 3/4 and 7/8 of the error vectors associated 
with the highest eigenvalues are reduced to about 1% of their initial values by 3 
nonstationary Point-Jacobi sweeps with 
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h — l/| A| maz , 1/ ^ (| A| maz ), 1 /( ^ I A | max) (8.2.1 ) 

if the eigenvalues are fairly evenly distributed between 0 and A maz . This leads to 
the concept of selectively annihilating clusters of eigenvectors from the error terms 
as a part of a total iteration process. This is also the basis for the success of the 
mixing and multigrid methods discussed below. 

8.3 Eigenvector and Eigenvalue Identification with Space Frequencies 

Consider the eigensystem of the model matrix |). The eigenvalues 

and eigenvectors for M — 5 are given in eqs (7.9.6) and (7.9.5), respectively. Notice 
that as the magnitudes of the eigenvalues increase the space frequency (number of 
sign changes) of the corresponding eigenvectors also increase. This has a rational 
explanation from the origin of the banded matrix. Note that 

(8.3.1) 


d 2 


dx 2 


- sin(mx) = — rrr sin(mx) 


and recall that 


X l 0 (8.3.2) 

We have seen that X~ x <j) represents a sine transform, and X<p, a sine synthesis. 
Therefore, the operation -~-zD(\) represents the numerical approximation of the 
multiplication of the appropriate sine wave by the negative square of its frequency, 
-m 2 . One finds that 


f>xx4> = “ 2 


X 


a>' ! 


1 

Ax 2 



-24 2 cos 


rmr 

M 4 1 


« — m 2 , rn << M 


(8.3.3) 


Hence, the correlation of large magnitudes of A m with high space-frequencies (see 
e.g., eq (7.6.1a)) is to be expected for these particular matrix operators. However, 
this correlation is not necessary in general. In fact, the complete counterexample of 
the above association is contained in the eigensystem for B( \ , 1, |). For this matrix 
one finds, from Section 5, exactly the opposite behavior. This is illustrated in eq 
(8.3.4) 
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X = 


1 

2 




2 


1 


iZI 

2 


1 

2 


iZi i iZI i 

2 1 2 2 

iZl o _ y 2* _ Vz 

2 u 2 2 

0-1 0 1 

_ n _ iZl 

2 u 2 2 

_ ^3 i _ ^3 1 


2 * 2 2 J 


( 8 . 3 . 4 ) 


r r — 

, corresponding 

I eigenvalue 


x l x 2 
1 . 866 ... 1.5 


*3 

1 


X 4 x 5 

0.5 . 134 ... 


Notice the eigenvectors with the lowest space-frequencies correspond to the eigen- 
values with maximum amplitude. Since the matrix structures that we are working 
with are subject to quite arbitrary preconditioning, this “inverse” structure can 
have practical implication. 
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9. EIGENSYSTEM MIXING 


The eigenvectors of 5(1,6, 1) are given by 


(l;m) = Sin 


3 


mn \ 

M + l) 


(9.1) 


and their structure does not depend upon the value of b. Therefore, the eigenvectors 
of the product 5(1,6], 1)5(1,62? 1) are also given by eq (8.2.1), and in general 


(Xjm) 



for 


k 


n*(iA,D 

k= 1 


(9.2) 


However, the eigenvalues of 5(1,6, 1) do depend upon 6, see eq (5.1.2), and when 
matrices have a common set of eigenvectors, the eigenvalues of their product are 
the product of their eigenvalues. Therefore, 


A»=nk-*~(i?TT)l/~n*<> ,6*,1) (9.3) 

k=\ 1 - ' + > k = 1 


This provides a way for constructing a process that modifies a matrix without 
changing its eigenvectors while changing the identities of the eigenvalues attached 
to them. Such a process is one way to mix the eigenvalue eigenvector structure 
of an iteration procedure so that a nonstationary method can remove more of the 
error without amplifying any of the vectors during the calculation. 

A simple example illustrates the point. Consider the linear averaging operator 
given by 


4> } - 4>j -i + 4>j + 4>j + 1 (9.4) 

when written in point operator form. If used twice on the 1-D model Dirichlet 
equation in a nonstationary Point Jacobi sequence, we have 


k + i =k + l, i, i ) 2 [b(i, - 2, i)^ n - f n 


(9.5) 


the A m eigenvalue spectrums of the matrices used above are given analytically by 
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(9.6) 


A m = — 2(1 — cos 6 m \ for 5(1, -2,1) 

A m = —2(1 - cos(30 m )j for 5(1, 1, 1) 2 £(1, - 2, 1) 
where 6 m = mir/(M -f 1) 


and are shown in Fig. 9.1. Remember that the eigenvectors of the two matrices 
indexed, in both cases, from 1 to M are identical. It is interesting to note that for 


M = 7 


5(l,l,l) 2 5(l,-2,l) = 


- 2-1 0 1 

- 1-2 0 0 1 
0 0-2 0 0 1 

1 0 0-2 0 0 1 

1 0 0-2 0 0 

1 0 0 - 2-1 
1 0 - 1-2 


(9.7) 


1 m M + 1 



Figure 9.1 - A m spectrums of 5(1, -2, 1) and 5(1, 1, l) 2 times B( 1, -2, 1) for 
odd M . 
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It is not the purpose of this report to invent optimum relaxation techniques, 
but rather to present some fundamental concepts by which they can be constructed. 
However, it is instructive to consider two ways in which the model equation 


B( l,-2,l)<£- / = 0 


(9.8) 


can be relaxed. First by a nonstationary Point-Jacobi or Richardson method given 
by 


6 


Ti -r i 



(9.9) 


which we refer to by PJ(h n ), and then by a combination of this method and the 
method given in eq (9.5) which we refer to by VM(h n ). 

Suppose the error content in the initial guess for <t> in eq (9.8) is expressed in 
terms of the amplitudes of the eigenvectors in (9.1), and suppose that initially each 
amplitude is given unit weight. That is 


M 

M x ) = X^ sin l m 0 Ax )l 

m— 1 

where x — j Ax and Ax = n/(M + 1) , and (9.10) 

2 M 

Mm) = <Mx) sin | j(mAz)) = (?) 

3 = 1 

Apply the Point-Jacobi sequence three times forming 

The error content vs. eigenvector number (in the sense defined in eq (7.6.1a)) is 
shown after each iteration in Fig. 9.2, where <j> represents ^ in vector space. This 
is, of course, very similar to Fig. 7.4.1 for Richardson's 3- step method using the 
Chebyshev h n sequence, except here the ordinate is the wave number. 
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Now consider a combination of the PJ and YM sequences given by 


= (9 - 12) 


where the initial conditions are again given by eq (9.10). For the first VM sweep 
with h - ] 4 the amplitudes of the eigenvectors associated with eigenvalues close 
to 4 are greatly reduced. These occur at m - M and M/3. These vicinities are 
further reduced by the second VM sweep with h = 1/3. Spikes of high amplitudes 
still remain at m = 0 and 2A//3. The spike at 2M/3 is easily eliminated by the two 
PJ sweeps using h = 1/3 since around m = 2M/3 the eigenvalues of the PJ matrix 
are -3. see Fig. 9.3. 
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The results of the two processes are compared in Fig. 9.4. A significant portion 
of the lower frequency content has been removed from the error content by the 
mixing strategy. The additional work required to do this would have to be compared 
with alternative methods for doing the same thing; but remember, the constraint 
imposed in constructing the results in Figures 9.2, 9.3, and 9.4 was that at no time 
was there any amplification of any eigenvector in the error content. 
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10. MULTIGRID STRATEGIES 


The idea of systematically using sets of coarser grids to accelerate the conver- 
gence of iteration schemes that arise from the numerical solution to partial differ- 
ential equations was made popular in this country by the work of Brandt, see ref. 7. 
There are many variations of the process, which is by no means unique, and many 
viewpoints of the underlying theory. The viewpoint presented here is a natural 
extension of the concepts discussed above. 

First of all we assume that the difference equations representing the basic 
partial differential equations are in a form that can be related to a matrix which 
has certain basic properties. This form can be arrived at “naturally” by simply 
replacing the derivatives in the PDE with difference schemes, as in the example 
given by eq (1.2.2), or it can be “contrived” by further conditioning, as in the 
examples given by eqs (2.2.9) and (2.2.10). These basic properties are: 


(1) The eigenvalues, A w , of the matrix are all real and negative. 

(2) The \ m are fairly evenly distributed between their maximum and 

minimum values. (10 1) 

(3) The eigenvectors associated with the eigenvalues having largest 
magnitudes can be correlated with high frequencies on the dif- 
ferencing mesh. 

These conditions are sufficient to ensure the validity of the process described next. 

Having preconditioned (if necessary) the basic finite differencing scheme by a 
procedure equivalent to the multiplication by a matrix C, we are led to the starting 
formulation 


C\A h 4 > oo - f b \ = 0 (10.2) 

where the matrix formed by the product CAb has the three properties in (10.1). In 
eq (10.2) the vector f b represents the boundary conditions and the forcing function, 
if any, and <7^ is a vector representing the desired exact solution. We start with 
some initial guess for 4 >oo and proceed through n iterations making use of some 
iterative process that greatly reduces the amplitudes of the eigenvectors associated 
with the eigenvalues in the range between |A| maz and 2|A| ma3; - We do not attempt 
to develop an optimum procedure here, but for clarity we suppose that the three- 
step Richardson method illustrated in Fig. 7.4.1 is used. At the end of the three 
steps we find r, the residual, where 
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r = C[A b ^-f b \ 


(10.3) 


We recall (Section 6) that the 4> used to compute r is composed of the exact solution 
< p qq and the error e in such a way that 


At - r = 0 

where A = CAt, 

If one could solve eq (10.4) for e then 


(10.4) 


**, = *- ' (10.5) 

Now' we can write the exact solution for e in terms of the eigenvectors of A, and the 
a eigenvalues of the Richardson process in the form (10.6) 


M/2 3 M3 

e = c m x m n "h ^ ^ 1 1 \o[X rn h m ^] (10.6) 

m- 1 n = l m=M/2 + 1 n=l 

' «*- 

very low amplitude 


Combining the properties of the Richardson algorithm and the conditions in (10.1), 
we can be sure that the high frequency content of e has been greatly reduced (about 
\% or less of its original value in the initial guess). 


Next we construct a permutation matrix which separates a vector into two 
parts, one containing the odd entries, and the other the even entries of the original 
matrix (or any other appropriate sorting which is consistent with the interpolation 
approximation to be discussed below). For example 


C2 


'o 

1 

0 

0 

0 

0 

o' 


ei 

e 4 


0 

0 

0 

1 

0 

0 

0 


C2 

«6 
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0 

0 

0 

0 

1 
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es 

Cl 

— 

1 

0 

0 

0 

0 

0 
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t A 

ez 
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0 

1 

0 

0 

0 

0 
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0 

0 

1 

0 

0 


ec 

_e 7 J 


0 

0 

0 

0 

0 

0 

1 


_e 7 _ 


e e 

eo 


= Pe 


(10.7) 
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Multiply eq (10.4) from the left by P and, since a permutation matrix has an inverse 
which is its transpose, we can write 


PA\P~ 1 P\e = Pr (10.8) 

The operation PAP~ l partitions the A matrix to form 

A\ A 2 
Az A 4 

Notice that 

A\e e + A 2 C 0 = r e (10.10) 


CO 


To 


(10.9) 


is an exact expression. At this point we make our one crucial assumption. It is 
that there is some connection between e e and e 0 brought about by the smoothing 
property of the Richardson relaxation procedure. Since the top half of the frequency 
spectrum has been removed, it is reasonable to suppose that the odd points are the 
average of the even points. For example 

e\ » l(e a + e 2 ) 

£3 ~ " ( e 2 + e 4) 

es « ^(e 4 + ee) 

e 7 ~ (®6 + e b ) 

It is important to notice that e a and ej, represent errors on the boundaries 
where the error is zero if the boundary conditions are given. It is also important to 
notice that we are dealing with the relation between e and r so the original boundary 
conditions and forcing function (which are contained in / in the basic formulation) 
no longer appear in the problem. Hence, no aliasing of these functions can occur in 
subsequent steps. Finally, notice that, in this formulation, the averaging of e is our 
only approximation, no operations on r are required or justified. 

If the boundary conditions are Dirichlet, e a and ef, are zero, and one can write 
for the example case 


or e 0 


A 2 e e 


( 10 . 11 ) 
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( 10 . 12 ) 


'1 0 O' 

. 1110 
* 2 ~ 2 0 1 1 

.0 0 1 . 

With this approximation eq (10.8) reduces to 


A c e e - r e = 0 

where A c ' [A i + ^ 2 ^ 2 ] 


(10.13) 


The form of A c , the matrix on the coarse mesh, is completely determined by the 
choice of the permutation matrix and the interpolation approximation. If the orig- 
inal A had been 5(1, -2, 1), our 7-point example would produce 


- 2 11 

2 1 1 
- 2 11 

11 - 2 
11 - 2 
11 - 2 
1 - 2 

and eq (10.13) gives 




(10.15) 

This process is deceptively simple. We started with the equation 5(1, -2, l)e = 
r on the fine mesh and reduced the problem to the equation ^5(1,— 2, l)e e — r e 
on the next coarser mesh. It appears as if the data on the odd points had been 
ignored altogether and a scaling factor had arbitrarily appeared. Such is not the 
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case, however, and except for the assumption in eq (10.11) the process is quite 
rigorous. 

If the boundary conditions are mixed Dirichlet Neumann, A in the 1-D model 
equation is 5(1,6, 1) where 6 = | —2, -2, ..., -2, — l] 7 . The eigensystem is given by 
eq (5.6.2). It is easy to show that the high space frequencies still correspond to the 
eigenvalues with high magnitudes, and, in fact, all of the conditions in (10.1) are 
met. However, the eigenvector structure is different from that shown in Fig. 7.6.1 
for Dirichlet conditions. In the present case they are given by 


Zjm 


= sin 



(2m - 1)tt V 

2M + 1 ) j ; 


m = 1,2, • • • , M 


(10.16a) 


and are illustrated in Fig. 10.1. All of them go through zero on the left (Dirichlet) 
side, and all of them reflect on the right (Neumann) side, being symmetrical about 
the point m -■ + ^ where x - n and their magnitude is 1. 



Figure 10.1 - Eigenvectors for the mixed Dirichlet-Neumann case. 


For Neumann conditions, the interpolation formula in eq (10.11) must be 
changed. In the particular case illustrated in Fig. 10.1, ej, is equal to e\j. If 
Neumann conditions are on the left, e a = ej. When et, = e m, the example in eq 
(10.12) changes to 
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(10.16b) 


1 0 O' 

1 1 0 
0 1 1 
0 0 2 , 

The permutation matrix remains the same and both A\ and A-i in the partitioned 
matrix PAP 1 are unchanged (only A 4 is modified by putting -1 in the lower right 
element). Therefore, we can construct the coarse matrix from 



A 1 


a 2 

/- 

A f 

N 

Ac 

- 2 

S 

1 1 

■v 

1 

'i 

1 1 


- 1 1/2 

- 2 

- 2 

+ 

1 1 
1 1 

’2 

1 1 
2. 


1/2 - 1 1/2 
1/2 - 1/2 


( 10 . 17 ) 

which gives us what we might have “expected” and shows us that the process is 
recursive. 

The remaining steps required to complete an entire multigrid process are rela- 
tively straightforward, but they vary depending on the problem and the user. The 
reduction can be, and usually is, carried to even coarser grids before returning to 
the finest level. However, in each case the appropriate permutation matrix and 
the interpolation approximation define both the down- and up-going paths. The 
details of finding optimum technique are, obviously, quite important but they are 
not discussed here. 
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